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The distribution function of suprathermal ions is found to be self-similar under conditions rel¬ 
evant to inertial confinement fusion hot-spots. By utilizing this feature, interference between the 
hydro-instabilities and kinetic effects is for the first time assessed quantitatively to find that the 
instabilities substantially aggravate the fusion reactivity reduction. The ion tail depletion is also 
shown to lower the experimentally inferred ion temperature, a novel kinetic effect that may ex¬ 
plain the discrepancy between the exploding pusher experiments and rad-hydro simulations and 
contribute to the observation that temperature inferred from DD reaction products is lower than 
from DT at National Ignition Facility. 


Recent exploding pusher experiments [lH6| reveal sub¬ 
stantial kinetic effects on the implosion performance. 
Specific mechanisms potentially responsible for these ob¬ 
servations include the inter-ion-species diffusion 7 10] 
and reactivity reduction due to ion tail depletion 11 -Il8| . 
Theoretical evaluation of these phenomena is challeng¬ 
ing, however, and while fully kinetic simulations allow 
study of a certain stage of implosion in specific config¬ 
urations 00, such calculations are computationally 
prohibitive for modeling of a realistic inertial confine¬ 
ment fusion (ICF) experiment. A substantial simplifi¬ 
cation results from treating thermal and suprathermal 
ions separately [H]. For the former, the mean-free-path 
A<°) is often much smaller than the characteristic scales 
of the system L, making fluid equations (including inter¬ 
species diffusion) a valid model. The latter constitute 
only a small fraction of the ion density, momentum, and 
energy and do not appear explicitly in the fluid equa¬ 
tions. However, it is the suprathermal ions that are most 
likely to undergo fusion reactions, so they do affect the 
fluid equations implicitly as an energy source. For these 
ions, the mean-free-path is much larger than A*- 0 ^ and 
can be comparable to L even if -C L. Hence, self- 
consistent modeling of ICF implosions would appear to 
require a kinetic treatment of suprathermal ions capable 
of predicting the fusion reactivity at each time step of 
the fluid equations’ evolution. 


While suprathermal ions can be described by a reduced 
linear (as opposite a to fully nonlinear) kinetic equa¬ 
tion [22| , this task is still non-trivial. All prior studies 
rely on either direct numerical solution |15Ml8l ] or phe¬ 
nomenological assumptions that affect the structure of 
the kinetic equation 3, 14]. Until now, no simple solu¬ 
tion to first-principles kinetic equation for the suprather¬ 
mal ions has been found even in the one-dimensional 
(ID) planar case. The issue becomes particularly press¬ 
ing in light of hydro-instabilities at the fuel-pusher inter¬ 
face |23H30(. It is near this interface that the suprather¬ 


mal ion distribution is modified most, so one should ex¬ 
pect substantial interference between the instabilities and 
the fusion reactivity. However, applying direct numeric 
methods to a complicated geometry is quite difficult and 
quantitative assessment of this interference has not been 
presented. 

In this Letter we demonstrate a physically intuitive, 
semi-analytical solution to the first-principles kinetic 
equation for suprathermal ions. This results from the 
self-similar structure of the ion distribution, scaling with 
the distance to the interface relative to the square of 
the ion energy. In the ID planar geometry, the solution 
agrees precisely with direct numerical results. Further¬ 
more, comparison with the numeric solution for the ID 
spherical geometry shows that the self-similar structure 
is robust against perturbations of the interface from the 
planar geometry. This allows us for the first time to 
evaluate the impact of hydro-instability on the reactivity 
reduction, which is found to be substantially enhanced. 
We also obtain a novel kinetic prediction for ICF exper¬ 
iments: that ion tail depletion results in the experimen¬ 
tally inferred temperature being lower than the actual 
one. 

We consider a spherically symmetric hot-spot with the 
radius Rh surrounded by a cold pusher. From symmetry, 
the distribution function f a of ion species a depends only 
on the three variables: radial coordinate r, particle speed 
v, and pitch angle 9 between velocity and radius vectors. 
Defining p, = cos 9 , we obtain the time-stationary Vlasov 
operator 
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The collision operator for species a is C„{/ a } = 
J2/3 C a p {fa}, where C a p denotes collisions of ion species 
a with ion species /3. We neglect ion-electron collisions 
and note that the suprathermal ions mostly collide with 
thermal ions, which are close to Maxwellian everywhere 
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outside a narrow vicinity of the boundary. Therefore, for 
suprathermal ions of species a, [22l. IffiD] 
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where m a , T a and ura = y2 Tajrria denote the particle 
mass, bulk temperature and thermal velocity of species 
a, respectively, and the collision frequency is defined by 


47 inpZlZje 4 

to>t« 


( 3 ) 


where Z a and n a are the charge number and the bulk ion 
density of species a, respectively, and In A is the Coulomb 
logarithm. 

The rate of energy exchange between thermal ions of 
species a and /3 is on the order of v a p ~ v aa ~ v 8p-, 
making their bulk temperatures equal, i.e. T a = To for 
all a. Assuming a flat temperature profile in the hot spot 
and equating the right sides of Eqs. © and (0 yield the 
stationary kinetic equation for the tail of f a 
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where x = r/Rh, £ = m a v 2 /2T$, 

N^ = VTa/iRhJftpVap), and G a = 

(Yftl3( m a/m/3)n/3Zjj / (j2p n pzf), with summa¬ 
tions over all ion species including a. The parameter 
G a depends only on the relative concentrations of the 
bulk ion densities; for example, in a 50/50 DT mixture, 
Gd = 5/6 and Gt = 5/4. The Knudsen number N j/ 1 ' 1 
is the ratio of the effective mean free path of a thermal 
ion aL 0) = VTa/{Yhf 3 v ap) and Rh and thus is the key 
parameter quantifying importance of the kinetic effects. 

The kinetic equation 0 needs to be accompanied by 
a condition at the interface between the hot and cold 
plasmas. A natural constraint results from assuming that 
there is no suprathermal ion inflow from the pusher into 
the hot-spot 


f a (x=l,-l<H<0,E) = 0. (5) 


Physically, one expects the solution to be effectively 
planar when the mean free path of a suprathermal ion 
with energy s, = £ 2 Xa \ is much less than Rh, or 

N^ ] e 2 < 1. (7) 

From direct numerical solution to the PDE problem 
formulated by Eq. 0 along with constraints © and 
© one can find that in this limit the distribution is 

self-similar, f a (x,/i,e) = f M <\>{ /"f 2 ,p), where f M = 

J\ K e 

n a (rn a /2tt T 0 f/ 2 e~ ~ e is Maxwellian. This finding has a 
transparent physical interpretation: for a given e devia¬ 
tion from equilibrium is controlled by the distance to the 
boundary y = Rh — r normalized to Aq^ . We then con¬ 
struct a solution that has this feature manifestly while 
keeping the pitch-angle scattering structure of the colli¬ 
sion operator 

fa = /M[l + ^c„^ n (p)e _<T " z ], (8) 


where 2 : = y/Xa^ = y/ (A^£ 2 ), c n are free constants, and 
eigenvalues and eigenvectors a n and il> n (/i) satisfy 
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Validity of this choice will be verified by comparing the 
resulting semi-analytical solution against the direct nu¬ 
merical one. We now proceed to identifying c n that are 
compatible with conditions © and ©. 

The spectrum of the operator of Eq. © is symmetric 
about zero (i.e. the eigenvalues come in pairs a n and 
—<J n ), a consequence of the right-left symmetry of the 
planar case. If the hot plasma occupies the half space 
—00 < y < 0, the constraint © dictates that only eigen¬ 
functions with <r n < 0 are included in the expansion ©. 
Imposing constraint © is a more non-trivial task, since 
it applies to /1 < 0 only. 

To implement this condition, we evaluate the matrix 
of the operator © over Legendre polynomials Tfc(/i). 
Solving the resulting eigenvalue problem gives if n {/i) = 
Ytk a nkPk{/f)- Defining b k = 1 + J2 n c n a nk for k = 0 
and bk = Y/ n c nd n k otherwise, condition © establishes 
a matrix relation between the vectors of odd and even 
coefficients bk 


&2fc+l = ^ D. 




( 10 ) 


In addition, the distribution function must be isotropic 
at the center due to symmetry 

df a {x = 0,n,e)/d(j, = 0. (6) 


where D nk = (4 k + 3)x 2 fc,2n+i with j32] 


(_l)( i +i+ 1 )/ 2 j!j! 
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Finally, inside the hot-spot we expect f a to become 
Maxwellian as £ approaches 1 from above since thermal 
ions are assumed to be close to equilibrium. 


Importantly, if N Legendre polynomials are kept in the 
expansion, this relation is equivalent to N/2 scalar equa¬ 
tions, which is exactly equal to the number of unknowns 
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in Eq. flS]) after eliminating modes with positive eigen¬ 
values [33|. Once the eigenvectors a n k are calculated, 
Eq. (gUD gives c n and therefore the distribution function 
through Eq. ((5J). 



FIG. 1. Spatial dependence of the pitch-angle averaged dis¬ 
tribution function relative to Maxwellian for several values 
of energy, as obtained from the direct numerical solution of 
Eq. ([4]) (solid) and the semi-analytical solution (dashed). 

To verify that this solution is precise in the planar case 
and investigate its robustness against deviations of the in¬ 
terface from the planar geometry, we compare it with the 
distribution function obtained by direct numerical solu¬ 
tion of Eq. (j4j. FigQshows a comparison of the deuteron 
distribution function relative to Maxwellian in the 50/50 
DT mixture with TV ^ = 0.01 where solid curves are to 
the numerical solution and dashed lines are to the planar 
solution (JH) with (1 — x)/(Nk£ 2 ) substituted for z. For 
energies such that Nk£ 2 <C 1, the solutions agree very 
well, which is consistent with formal planar limit condi¬ 
tion Q. Furthermore, the agreement is good even for 
Nk£ 2 ~ 0.5 (red curve). Comparisons for various TVj^ 
and mixture compositions show similar agreement. 

When A o' 1 = £ 2 Aq 0 ^ > Rh, the planar solution over¬ 
estimates the suprathermal ion population, reflecting a 
shortcoming of the planar model in the context of spher¬ 
ical geometry when ions “feel” not only the distance to 
the boundary, but also the distance to the center of the 
hot-spot. To investigate significance of this spherical ef¬ 
fect we consider the reactivity (av). 



FIG. 2. Volume averaged DT reactivity relative to 
Maxwellian computed from the direct numerical solution of 
Eq. (solid) and from the semi-analytical solution (dashed). 


The fusion cross-section a is defined in absolute en¬ 
ergy units, while f a deviation from Maxwellian (as well 
as deviation between planar and spherical solutions) is 
governed by the energy £ normalized to the main ion 
temperature. Thus, for a fixed TV ^ and higher T 0 the 
reactivity is due to lower e, thereby diminishing these 
deviations. Conversely, for fixed To, larger corre¬ 

sponds to a larger mean-free-path for any given e. making 
the effect stronger. This qualitative picture is supported 
by the results in Fig [2j which presents the volume aver¬ 
aged reactivity (av) = V -1 f v d 3 r (av) for the DT fusion 
reaction in the 50/50 DT mixture. For T 0 > 5 keV, the 
predictions of the planar model and the direct numeri¬ 
cal calculation are within 15% even for TV ^ = 0.1 and 
within 5% for TV ^ < 0.03. In the 1-5 keV range the dis¬ 
crepancy is larger, but still within 12% for TVj^ = 0.01 
and 0.03, for which deviation from ID case will be con¬ 
sidered. Hence, with respect to reactivity, the distance 
to the boundary is the only relevant scale and the re¬ 
sulting planar solution is robust to perturbations of the 
boundary from the planar geometry. 


The above allows us for the first time to perform a 
quantitative assessment of kinetic effects on the reactiv¬ 
ity in the presence of hydro-instability. The entire implo¬ 
sion analysis would include modeling the instability along 
with the kinetic effects. For the purpose of demonstra¬ 
tion, here we consider a sample instability resulting in 
the interface perturbed as R = Rh + A R cos (to$) , where 
R and d are the radial and polar angle coordinates of the 
hot-spot boundary and symmetry in the azimuthal an¬ 
gle ip is assumed. Experiments and simulations indicate 
that A R/R can be as large as 1/3 HI HI El 0 and 
to imitate the spatial structure of the unstable region we 
take m = 20. Then, we evaluate the overall reactivity 
reduction as follows: for a point A inside the hot-spot 
the point B on the perturbed surface is found such that 
the distance AB is the shortest; the distribution function 
is evaluated from our solution with AB inserted for y\ lo¬ 
cal reactivity reduction at each point A is calculated and 
averaged over the perturbed volume V. This procedure 
is computationally inexpensive and can be applied read¬ 
ily to an arbitrarily perturbed surface. Fig [3] presents 
results of this calculation for the DT reaction in a 50/50 
DT mixture for A R/Rh equal to 1/5 (dashed line) and 
1/3 (dotted line) and Knudsen numbers of 0.01 and 0.03. 
Introduction of a perturbation even at the conservative 
level of A R/Rh = 1/5 can double the reactivity reduc¬ 
tion in the 1-10 keV range. Importantly, the change in 
reactivity from the surface perturbation is substantially 
larger than discrepancy between predictions of the nu¬ 
merical and semi-analytical solutions in Fig [2] for any 
given T 0 . 

Finally, we demonstrate a novel effect of the ion tail de¬ 
pletion on the experimentally inferred ion temperature. 
In ICF experiments the hot-spot ion temperature is de- 
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FIG. 3. Volume averaged DT reactivity relative to 
Maxwellian computed for the unperturbed spherical hot¬ 
spot boundary (solid) and for the boundary perturbed with 
A R/Rh = 1/5 (dashed) and A R/Rh = 1/3 (dotted). 


FIG. 4. Volume averaged DT and DD burn temperature rel¬ 
ative to the bulk ion temperature To computed from the di¬ 
rect numerical solution for N^ >) — 0.05 (solid), N^ = 0.1 
(dashed) and = 0.2 (dotted). 


duced from the width of the reaction products’ spectra, 
which can be related to the mean square of the center-of- 
mass velocity of ion pairs undergoing this reaction [341 


T — 

exp — 


TOi +TO2 / d 3 Vid 3 V 2 fif 2 (TVU\ 


c.m. 


f d 3 V 1 d 3 V 2 f 1 f 2 <JV 


( 11 ) 


where v = |u| = |ui — v 2 \ and u c . m . = (m\V\ + 
Tn 2 v 2 )/(mi + m 2 ). Since the reaction cross section a 
depends only on v, the expression on the right side 
of Eq. (fill can be shown to recover To when f\ }2 are 
Maxwellian regardless of the fusion reaction. For distri¬ 
bution functions different from /m, one should expect 
T exp to also be different from To and, equally impor¬ 
tantly, to depend on a. The DD cross-section is due 
to larger energies than the DT cross-section and, in the 
50/50 DT mixture, fjj is farther from equilibrium than 
fr due to Go < Gt ■ As a result, T exp associated with 
DD reaction is lower than that associated with DT reac¬ 
tion, which is confirmed by Fig [4] showing the two tem¬ 
peratures evaluated from the direct numerical solution 
according to Eq. m We also see that the reduction 
in T exp is less than that in (crv). In addition, it can be 
found that, unlike (crv) case, it is the spherical effects on 
the distribution function rather than its planar structure 
that reflect stronger on T exp . This difference between 
the (crv) and T exp , to leading order, results from the for¬ 
mer being governed by the ion number density within 
the Gamow window, whereas the latter is governed by 
higher moments and finer features of the distribution 
function. The size of the predicted effect is consistent 
with the discrepancy between the temperature observed 
in exploding pusher experiments, in which Nk can be 
even higher than shown in Fig|4[ and that predicted by 
simulations US]- Of course, in these experiments kinetic 
effects can also lower the actual temperature by reducing 
the shock heating; however, employing standard formu¬ 
lae for the spectrum width would diminish the inferred 
temperature even further. 

The DD temperature lower than DT by about 25% is 
often observed in cryo-implosions at NIF [35J. The only 


explanation available to date is based on the bulk fluid 
motion [36| and a detailed analysis of this effect in ampli¬ 
cation to turbulent implosions is presented in Ref. 1371 ]. 
However, to give the 25% difference between the burn 
temperatures this mechanism requires all the plasma en¬ 
ergy to be kinetic that is hardly possible during ICF burn. 
In a more realistic scenario with the turbulent energy 
on the order of the thermal energy it can give 10-15% 
only [13] • The following estimate demonstrates that the 
newly predicted effect can be responsible for the remain¬ 
ing part of the difference. 

The turbulence strongly distorts the mixing layer 
topology, so the approach employed for evaluating the 
role of instabilities no longer applies. To estimate the 
above kinetic effects in this situation one can view the 
mix layer as a suspension of hot plasma droplets in a 
cold pusher material [16]. The effective Knudsen num¬ 


ber associated with a given droplet is Xa 1 / L, where L 
is the characteristic radius of the droplet rather than of 
the hot-spot. The nominal Knudsen number for NIF, 
Xa' 1 /Rh, is a few percent 0 and Rh ^ T, so the effec¬ 
tive Knudsen number should be taken to be a fraction of 
unity. Results of Fig 0 ] give that the difference between 
the DD and DT temperatures grows with Nk, becom¬ 
ing about 12% in the 1-5 keV range for Nk = 0.2. For 
higher Knudsen numbers, which are likely according to 
the above, our approach to suprathermal ions is com¬ 
plicated; yet, the trend suggests that the DD and DT 
temperatures fall apart even further. Also, the mix layer 
occupies the outer radii and therefore constitutes a large 
fraction of the hot-spot. Hence, 10-15% of the difference 
between the DD and DT temperatures can be realisti¬ 
cally attributed to the ion tail depletion and the newly 
presented mechanism together with the bulk fluid mo¬ 
tion effect considered earlier [37| can fully explain the 
observations. 

Our predictions are limited by shortcomings of the ki¬ 
netic equation 0, which assumes a constant ion tem¬ 
perature in the hot-spot. To investigate robustness of 
the newly found effects with respect to a more realis- 
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tic temperature profile, we have conducted direct nu¬ 
merical simulations for an isobaric hot-spot with the 
main ion temperatu re g iven by commonly observed T = 
T 0 [ 1 — (r/fi/j) 2 ] 2 / 7 [39]. The results for the reactivity 
reduction and the experimentally inferred ion temper¬ 
atures in terms of the volume averaged Knudsen num¬ 
ber NjP' 1 and bulk ion temperature T turn out to be 
in reasonable agreement with the results of Figs [2] and 
|]|] and show the same general trends with these parame¬ 
ters. Furthermore, while the flat profile model presented 
in this Letter slightly overpredicts the reactivity reduc¬ 
tion, it underpredicts the difference between the exper¬ 
imentally inferred DD and DT temperatures. We thus 
expect that our conclusions qualitatively persist in prac¬ 
tical hot-spot configurations, though some quantitative 
changes are possible. 

To summarize, a semi-analytical solution to first- 
principles equation for the suprathermal ions in ID ge¬ 
ometry has been constructed and provides a computa¬ 
tionally expedient tool for investigating kinetic effects 
in complicated geometries. In particular, the analy¬ 
sis demonstrates that hydrodynamic instabilities at hot¬ 
spot/pusher interfaces can substantially aggravate the 
reactivity reduction. Moreover, the ion tail depletion re¬ 
sults in the experimentally inferred core ion temperatures 
being lower than the actual ones, which may explain re¬ 
cent measurements in exploding pusher implosions and 
contribute to the observation that DD burn temperature 
is lower than DT burn temperature at NIF. 
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